Libraries
library(survival)
library(FRESA.CAD)
## Loading required package: Rcpp
## Loading required package: stringr
## Loading required package: miscTools
## Loading required package: Hmisc
##
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:base':
##
## format.pval, units
## Loading required package: pROC
## Type 'citation("pROC")' for a citation.
##
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
##
## cov, smooth, var
#library(corrplot)
#source("~/GitHub/FRESA.CAD/R/RRPlot.R")
#source("~/GitHub/FRESA.CAD/R/PoissonEventRiskCalibration.R")
op <- par(no.readonly = TRUE)
pander::panderOptions('digits', 3)
#pander::panderOptions('table.split.table', 400)
pander::panderOptions('keep.trailing.zeros',TRUE)
The data set
data(cancer)
colon <- subset(colon,etype==1)
colon$etype <- NULL
rownames(colon) <- colon$id
colon$id <- NULL
colon <- colon[complete.cases(colon),]
time <- colon$time
status <- colon$status
data <- colon
data$time <- NULL
data$study <- NULL
table(data$status)
0 1 442 446
dataColon <- as.data.frame(model.matrix(status~.*age,data))
dataColon$`(Intercept)` <- NULL
dataColon$time <- time/365
dataColon$status <- status
colnames(dataColon) <-str_replace_all(colnames(dataColon),":","_")
colnames(dataColon) <-str_replace_all(colnames(dataColon),"\\.","_")
colnames(dataColon) <-str_replace_all(colnames(dataColon),"\\+","_")
data <- NULL
trainsamples <- sample(nrow(dataColon),0.7*nrow(dataColon))
dataColonTrain <- dataColon[trainsamples,]
dataColonTest <- dataColon[-trainsamples,]
pander::pander(table(dataColonTrain$status))
pander::pander(table(dataColonTest$status))
Cox Model Performance
Here we evaluate the model using the RRPlot() function.
The evaluation of the raw Cox model with RRPlot()
Here we will use the predicted event probability assuming a baseline
hazard for events withing 5 years
index <- predict(ml,dataColonTrain)
timeinterval <- 2*mean(subset(dataColonTrain,status==1)$time)
h0 <- sum(dataColonTrain$status & dataColonTrain$time <= timeinterval)
h0 <- h0/sum((dataColonTrain$time > timeinterval) | (dataColonTrain$status==1))
rdata <- cbind(dataColonTrain$status,ppoisGzero(index,h0))
rrAnalysisTrain <- RRPlot(rdata,atRate=c(0.90),
timetoEvent=dataColonTrain$time,
title="Raw Train: Colon Cancer",
ysurvlim=c(0.00,1.0),
riskTimeInterval=timeinterval)






Uncalibrated Performance Report
pander::pander(t(rrAnalysisTrain$keyPoints),caption="Threshold values")
Threshold values
| Thr |
0.4610 |
0.344 |
0.292 |
0.27013 |
0.4997 |
| RR |
1.6114 |
1.623 |
2.170 |
1.00000 |
1.5921 |
| RR_LCI |
1.3949 |
1.389 |
1.319 |
0.00000 |
1.3617 |
| RR_UCI |
1.8615 |
1.897 |
3.570 |
0.00000 |
1.8615 |
| SEN |
0.2785 |
0.560 |
0.962 |
1.00000 |
0.1772 |
| SPE |
0.8951 |
0.685 |
0.121 |
0.00656 |
0.9410 |
| BACC |
0.5868 |
0.623 |
0.542 |
0.50328 |
0.5591 |
| NetBenefit |
0.0976 |
0.204 |
0.312 |
0.32829 |
0.0612 |
pander::pander(t(rrAnalysisTrain$OERatio$estimate),caption="O/E Ratio")
O/E Ratio
| 0.98 |
0.875 |
1.09 |
0.759 |
pander::pander(t(rrAnalysisTrain$OE95ci),caption="O/E Mean")
O/E Mean
| 1.48 |
1.48 |
1.46 |
1.51 |
pander::pander(t(rrAnalysisTrain$OAcum95ci),caption="O/Acum Mean")
O/Acum Mean
| 1.39 |
1.39 |
1.39 |
1.4 |
pander::pander(rrAnalysisTrain$c.index$cstatCI,caption="C. Index")
pander::pander(t(rrAnalysisTrain$ROCAnalysis$aucs),caption="ROC AUC")
ROC AUC
| 0.663 |
0.621 |
0.706 |
pander::pander((rrAnalysisTrain$ROCAnalysis$sensitivity),caption="Sensitivity")
Sensitivity
| 0.278 |
0.23 |
0.331 |
pander::pander((rrAnalysisTrain$ROCAnalysis$specificity),caption="Specificity")
Specificity
| 0.898 |
0.859 |
0.93 |
pander::pander(t(rrAnalysisTrain$thr_atP),caption="Probability Thresholds")
Probability Thresholds
| 0.462 |
pander::pander(t(rrAnalysisTrain$RR_atP),caption="Risk Ratio")
Quitting from lines 109-122 (Colon.Rmd) Error in
t.default(rrAnalysisTrain$RR_atP) : argument is not a matrix Calls:
… eval_with_user_handlers -> eval -> eval ->
-> t -> t.default In addition: Warning messages: 1: In
doTryCatch(return(expr), name, parentenv, handler) : “quiet” is not a
graphical parameter 2: In doTryCatch(return(expr), name, parentenv,
handler) : “quiet” is not a graphical parameter 3: In
doTryCatch(return(expr), name, parentenv, handler) : “quiet” is not a
graphical parameter 4: In doTryCatch(return(expr), name, parentenv,
handler) : “quiet” is not a graphical parameter 5: In
doTryCatch(return(expr), name, parentenv, handler) : “quiet” is not a
graphical parameter 6: In doTryCatch(return(expr), name, parentenv,
handler) : “quiet” is not a graphical parameter 7: In
doTryCatch(return(expr), name, parentenv, handler) : “quiet” is not a
graphical parameter 8: In doTryCatch(return(expr), name, parentenv,
handler) : “quiet” is not a graphical parameter 9: In
doTryCatch(return(expr), name, parentenv, handler) : “quiet” is not a
graphical parameter
pander::pander(rrAnalysisTrain$surdif,caption="Logrank test")
Logrank test Chisq = 58.848822 on 1 degrees of freedom, p =
0.000000
| class=0 |
502 |
228 |
274.1 |
7.74 |
58.8 |
| class=1 |
119 |
88 |
41.9 |
50.60 |
58.8 |
Cox Calibration
op <- par(no.readonly = TRUE)
calprob <- CoxRiskCalibration(ml,dataColonTrain,"status","time")
pander::pander(c(h0=calprob$h0,
Gain=calprob$hazardGain,
DeltaTime=calprob$timeInterval),
caption="Cox Calibration Parameters")
Calibrated Train Performance
pander::pander(t(rrAnalysisTrain$keyPoints),caption="Threshold values")
Threshold values
| Thr |
0.610 |
0.474 |
0.409 |
0.38132 |
0.500 |
| RR |
1.611 |
1.623 |
2.122 |
1.00000 |
1.544 |
| RR_LCI |
1.395 |
1.389 |
1.292 |
0.00000 |
1.334 |
| RR_UCI |
1.862 |
1.897 |
3.485 |
0.00000 |
1.787 |
| SEN |
0.278 |
0.560 |
0.962 |
1.00000 |
0.418 |
| SPE |
0.895 |
0.685 |
0.118 |
0.00656 |
0.787 |
| BACC |
0.587 |
0.623 |
0.540 |
0.50328 |
0.602 |
| NetBenefit |
0.061 |
0.146 |
0.190 |
0.20817 |
0.108 |
pander::pander(t(rrAnalysisTrain$OERatio$estimate),caption="O/E Ratio")
O/E Ratio
| 0.779 |
0.695 |
0.87 |
4.3e-06 |
pander::pander(t(rrAnalysisTrain$OE95ci),caption="O/E Mean")
O/E Mean
| 0.988 |
0.988 |
0.974 |
1 |
pander::pander(t(rrAnalysisTrain$OAcum95ci),caption="O/Acum Mean")
O/Acum Mean
| 1.06 |
1.06 |
1.05 |
1.06 |
pander::pander(rrAnalysisTrain$c.index$cstatCI,caption="C. Index")
pander::pander(t(rrAnalysisTrain$ROCAnalysis$aucs),caption="ROC AUC")
ROC AUC
| 0.663 |
0.621 |
0.706 |
pander::pander((rrAnalysisTrain$ROCAnalysis$sensitivity),caption="Sensitivity")
Sensitivity
| 0.278 |
0.23 |
0.331 |
pander::pander((rrAnalysisTrain$ROCAnalysis$specificity),caption="Specificity")
Specificity
| 0.898 |
0.859 |
0.93 |
pander::pander(t(rrAnalysisTrain$thr_atP),caption="Probability Thresholds")
Probability Thresholds
| 0.611 |
pander::pander(t(rrAnalysisTrain$RR_atP),caption="Risk Ratio")
Quitting from lines 161-174 (Colon.Rmd) Error in
t.default(rrAnalysisTrain$RR_atP) : argument is not a matrix Calls:
… eval_with_user_handlers -> eval -> eval ->
-> t -> t.default In addition: Warning messages: 1: In
doTryCatch(return(expr), name, parentenv, handler) : “quiet” is not a
graphical parameter 2: In doTryCatch(return(expr), name, parentenv,
handler) : “quiet” is not a graphical parameter 3: In
doTryCatch(return(expr), name, parentenv, handler) : “quiet” is not a
graphical parameter 4: In doTryCatch(return(expr), name, parentenv,
handler) : “quiet” is not a graphical parameter 5: In
doTryCatch(return(expr), name, parentenv, handler) : “quiet” is not a
graphical parameter 6: In doTryCatch(return(expr), name, parentenv,
handler) : “quiet” is not a graphical parameter 7: In
doTryCatch(return(expr), name, parentenv, handler) : “quiet” is not a
graphical parameter 8: In doTryCatch(return(expr), name, parentenv,
handler) : “quiet” is not a graphical parameter 9: In
doTryCatch(return(expr), name, parentenv, handler) : “quiet” is not a
graphical parameter
pander::pander(rrAnalysisTrain$surdif,caption="Logrank test")
Logrank test Chisq = 58.848822 on 1 degrees of freedom, p =
0.000000
| class=0 |
502 |
228 |
274.1 |
7.74 |
58.8 |
| class=1 |
119 |
88 |
41.9 |
50.60 |
58.8 |
Test Performance
pander::pander(t(rrAnalysisTest$keyPoints),caption="Threshold values")
Threshold values
| Thr |
0.6083 |
0.451 |
0.407 |
3.93e-01 |
0.501 |
| RR |
1.7584 |
3.608 |
8.772 |
1.98e+01 |
1.875 |
| RR_LCI |
1.4078 |
2.212 |
1.305 |
4.34e-02 |
1.486 |
| RR_UCI |
2.1962 |
5.886 |
58.958 |
9.01e+03 |
2.367 |
| SEN |
0.3231 |
0.892 |
0.992 |
1.00e+00 |
0.492 |
| SPE |
0.8905 |
0.489 |
0.117 |
2.92e-02 |
0.803 |
| BACC |
0.6068 |
0.691 |
0.555 |
5.15e-01 |
0.648 |
| NetBenefit |
0.0701 |
0.219 |
0.172 |
1.65e-01 |
0.138 |
pander::pander(t(rrAnalysisTest$OERatio$estimate),caption="O/E Ratio")
O/E Ratio
| 0.708 |
0.591 |
0.84 |
3.52e-05 |
pander::pander(t(rrAnalysisTest$OE95ci),caption="O/E Mean")
O/E Mean
| 0.931 |
0.932 |
0.907 |
0.956 |
pander::pander(t(rrAnalysisTest$OAcum95ci),caption="O/Acum Mean")
O/Acum Mean
| 1.07 |
1.07 |
1.07 |
1.08 |
pander::pander(rrAnalysisTest$c.index$cstatCI,caption="C. Index")
pander::pander(t(rrAnalysisTest$ROCAnalysis$aucs),caption="ROC AUC")
ROC AUC
| 0.741 |
0.682 |
0.799 |
pander::pander((rrAnalysisTest$ROCAnalysis$sensitivity),caption="Sensitivity")
Sensitivity
| 0.323 |
0.244 |
0.411 |
pander::pander((rrAnalysisTest$ROCAnalysis$specificity),caption="Specificity")
Specificity
| 0.898 |
0.834 |
0.943 |
pander::pander(t(rrAnalysisTest$thr_atP),caption="Probability Thresholds")
Probability Thresholds
| 0.611 |
pander::pander(t(rrAnalysisTest$RR_atP),caption="Risk Ratio")
Quitting from lines 196-209 (Colon.Rmd) Error in
t.default(rrAnalysisTest$RR_atP) : argument is not a matrix Calls:
… eval_with_user_handlers -> eval -> eval ->
-> t -> t.default In addition: Warning messages: 1: In
doTryCatch(return(expr), name, parentenv, handler) : “quiet” is not a
graphical parameter 2: In doTryCatch(return(expr), name, parentenv,
handler) : “quiet” is not a graphical parameter 3: In
doTryCatch(return(expr), name, parentenv, handler) : “quiet” is not a
graphical parameter 4: In doTryCatch(return(expr), name, parentenv,
handler) : “quiet” is not a graphical parameter 5: In
doTryCatch(return(expr), name, parentenv, handler) : “quiet” is not a
graphical parameter 6: In doTryCatch(return(expr), name, parentenv,
handler) : “quiet” is not a graphical parameter 7: In
doTryCatch(return(expr), name, parentenv, handler) : “quiet” is not a
graphical parameter 8: In doTryCatch(return(expr), name, parentenv,
handler) : “quiet” is not a graphical parameter 9: In
doTryCatch(return(expr), name, parentenv, handler) : “quiet” is not a
graphical parameter
pander::pander(rrAnalysisTest$surdif,caption="Logrank test")
Logrank test Chisq = 31.486584 on 1 degrees of freedom, p =
0.000000
| class=0 |
211 |
88 |
110.7 |
4.64 |
31.5 |
| class=1 |
56 |
42 |
19.3 |
26.54 |
31.5 |
Cross-Validation
Here we will cross validate the training set and evaluate also on the
testing set. The h0 and the timeinterval are the ones estimated on the
calibration process
rcv <- randomCV(theData=dataColonTrain,
theOutcome = Surv(time,status)~1,
fittingFunction=BSWiMS.model,
trainFraction = 0.75,
repetitions=50,
classSamplingType = "Pro",
testingSet=dataColonTest
)
.[++++].[++++].[+++-].[++++].[++++].[++++–].[++++-].[++++—].[++++–].[++++]10
Tested: 854 Avg. Selected: 6.4 Min Tests: 1 Max Tests: 10 Mean Tests:
4.953162 . MAD: 0.4796788
.[+++].[+++++].[++++-].[+++++].[++++-].[++++].[+++].[++++].[+++++].[++++]20
Tested: 887 Avg. Selected: 6.8 Min Tests: 1 Max Tests: 20 Mean Tests:
9.537768 . MAD: 0.4792116
.[++++].[+++–].[++++].[++++].[++++].[++++].[++–].[++++].[++++–].[+++–]30
Tested: 888 Avg. Selected: 6.666667 Min Tests: 1 Max Tests: 30 Mean
Tests: 14.29054 . MAD: 0.4791815
.[++–].[+++–].[++++].[+++++].[++++—].[++++++].[+++].[++++].[+++-].[++++-]40
Tested: 888 Avg. Selected: 6.85 Min Tests: 1 Max Tests: 40 Mean Tests:
19.05405 . MAD: 0.4789083
.[++++-].[+++++].[+++-].[++++].[+++-+-].[+++].[++-].[+++-].[++++].[++-]50
Tested: 888 Avg. Selected: 6.72 Min Tests: 3 Max Tests: 50 Mean Tests:
23.81757 . MAD: 0.4790776
stp <- rcv$survTestPredictions
stp <- stp[!is.na(stp[,4]),]
bbx <- boxplot(unlist(stp[,1])~rownames(stp),plot=FALSE)
times <- bbx$stats[3,]
status <- boxplot(unlist(stp[,2])~rownames(stp),plot=FALSE)$stats[3,]
prob <- ppoisGzero(boxplot(unlist(stp[,4])~rownames(stp),plot=FALSE)$stats[3,],h0)
rdatacv <- cbind(status,prob)
rownames(rdatacv) <- bbx$names
names(times) <- bbx$names
rrAnalysisCVTest <- RRPlot(rdatacv,atThr = rrAnalysisTrain$thr_atP,
timetoEvent=times,
title="CV Test: Colon Cancer",
ysurvlim=c(0.00,1.0),
riskTimeInterval=timeinterval)






CV Test Performance
pander::pander(t(rrAnalysisCVTest$keyPoints),caption="Threshold values")
Threshold values
| Thr |
0.6109 |
0.483 |
0.425 |
0.369 |
0.500 |
| RR |
1.5886 |
1.677 |
2.446 |
1.000 |
1.604 |
| RR_LCI |
1.4053 |
1.473 |
1.681 |
0.000 |
1.417 |
| RR_UCI |
1.7959 |
1.909 |
3.557 |
0.000 |
1.815 |
| SEN |
0.3229 |
0.538 |
0.951 |
1.000 |
0.435 |
| SPE |
0.8620 |
0.719 |
0.176 |
0.000 |
0.787 |
| BACC |
0.5924 |
0.629 |
0.564 |
0.500 |
0.611 |
| NetBenefit |
0.0543 |
0.140 |
0.174 |
0.211 |
0.113 |
pander::pander(t(rrAnalysisCVTest$OERatio$estimate),caption="O/E Ratio")
O/E Ratio
| 0.751 |
0.683 |
0.824 |
2.86e-10 |
pander::pander(t(rrAnalysisCVTest$OE95ci),caption="O/E Mean")
O/E Mean
| 0.936 |
0.936 |
0.925 |
0.947 |
pander::pander(t(rrAnalysisCVTest$OAcum95ci),caption="O/Acum Mean")
O/Acum Mean
| 1.02 |
1.02 |
1.02 |
1.02 |
pander::pander(rrAnalysisCVTest$c.index$cstatCI,caption="C. Index")
pander::pander(t(rrAnalysisCVTest$ROCAnalysis$aucs),caption="ROC AUC")
ROC AUC
| 0.668 |
0.633 |
0.704 |
pander::pander((rrAnalysisCVTest$ROCAnalysis$sensitivity),caption="Sensitivity")
Sensitivity
| 0.321 |
0.278 |
0.366 |
pander::pander((rrAnalysisCVTest$ROCAnalysis$specificity),caption="Specificity")
Specificity
| 0.862 |
0.826 |
0.893 |
pander::pander(t(rrAnalysisCVTest$thr_atP),caption="Probability Thresholds")
Probability Thresholds
| 0.611 |
pander::pander(t(rrAnalysisCVTest$RR_atP),caption="Risk Ratio")
Quitting from lines 251-263 (Colon.Rmd) Error in
t.default(rrAnalysisCVTest$RR_atP) : argument is not a matrix Calls:
… eval_with_user_handlers -> eval -> eval ->
-> t -> t.default In addition: Warning messages: 1: In
doTryCatch(return(expr), name, parentenv, handler) : “quiet” is not a
graphical parameter 2: In doTryCatch(return(expr), name, parentenv,
handler) : “quiet” is not a graphical parameter 3: In
doTryCatch(return(expr), name, parentenv, handler) : “quiet” is not a
graphical parameter 4: In doTryCatch(return(expr), name, parentenv,
handler) : “quiet” is not a graphical parameter 5: In
doTryCatch(return(expr), name, parentenv, handler) : “quiet” is not a
graphical parameter 6: In doTryCatch(return(expr), name, parentenv,
handler) : “quiet” is not a graphical parameter 7: In
doTryCatch(return(expr), name, parentenv, handler) : “quiet” is not a
graphical parameter 8: In doTryCatch(return(expr), name, parentenv,
handler) : “quiet” is not a graphical parameter 9: In
doTryCatch(return(expr), name, parentenv, handler) : “quiet” is not a
graphical parameter
pander::pander(rrAnalysisCVTest$surdif,caption="Logrank test")
Logrank test Chisq = 72.028528 on 1 degrees of freedom, p =
0.000000
| class=0 |
684 |
303 |
370.1 |
12.2 |
72 |
| class=1 |
204 |
143 |
75.9 |
59.3 |
72 |